收藏|limma转录组定量中文操作指南
很多人可能以为RNA-seq数据分析很难学，然而实际上并不是这样的。在limma-voom的作者Charity Law和她的同事们写的RNA-seq数据差异表达分析流程的教程中，RNA-seq数据分析被形容为“如同数1-2-3一样简单”。这篇教程中的所有分析都是在R语言中进行的，而且用到的所有包都是完全免费的。不仅如此，这篇教程最近由作者实验室的中国学生（也就是笔者本人2333）翻译成了中文 ，更为广大同胞免去了看英文教程时难以理解的痛苦！上网搜索“RNAseq123”并进入其Bioconductor包页面，即可查看教程中英文原文和其中所有步骤的代码。
这篇推文里，我们缩减了一下教程原文，简单地介绍了RNA-seq差异表达分析的数据整合、数据预处理、差异表达分析、热图和MD图的绘制以及camera基因集检验。如果想要进一步学习，更推荐大家查看教程的原文（戳“阅读原文”即可访问）。如果在学习的过程中遇到任何问题，都可以在Bioconductor support论坛提问。
1、背景介绍
这篇教程中用到的数据（GEO序列号: GSE63310）由雌性处女小鼠的三种乳腺细胞类型组成，即基底（basal）、管腔祖细胞（liminal progenitor, LP）和成熟管腔（mature luminal, ML），每组各设三个生物学重复。
文中所描述的数据分析是从已经经过比对和基因水平计数的表达计数（counts）数据起步的。如果你手头想要分析的数据是原始的序列也不要慌，作者推荐使用Rsubread包的align函数和featureCounts函数在R中进行比对和计数，这两个函数的用法也不难学。
首先载入分析需要用到的软件包。
library(limma)library(Glimma)
library(edgeR)
library(Mus.musculus)
2、数据整合
使用edgeR的readDGE函数读入数据，合并并创建DGEList对象x，其中包含一个计数矩阵，它的27179行分别对应每个基因的Entrez基因ID，九列分别对应此实验中的每个样品。如果来自所有样品的计数存储在同一个文件中，那么数据可以首先读入R再使用DGEList函数转换为一个DGEList对象。
"GSM1545538_purep53.txt", "GSM1545539_JMS8-2.txt",
"GSM1545540_JMS8-3.txt", "GSM1545541_JMS8-4.txt",
"GSM1545542_JMS8-5.txt", "GSM1545544_JMS9-P7c.txt",
"GSM1545545_JMS9-P8c.txt")
x <- readDGE(files, columns=c(1,3))
组织样品信息，将样品分组、处理、批次等与实验设计有关的样品信息与计数矩阵的列关联。
samplenames <- substring(colnames(x), 12, nchar(colnames(x)))colnames(x) <- samplenames
group <- as.factor(c("LP", "ML", "Basal", "Basal", "ML", "LP",
"Basal", "ML", "LP"))
x$samples$group <- group
lane <- as.factor(rep(c("L004","L006","L008"), c(3,4,2)))
x$samples$lane <- lane
使用Mus.musculus包，利用数据集中的Entrez基因ID来检索相对应的基因符号和染色体信息，将这些基因注释信息存储到DGEList对象中的第二个数据框genes。
genes <- select(Mus.musculus, keys=geneid, columns=c("SYMBOL", "TXCHROM"),
keytype="ENTREZID")
为了处理重复的基因ID，选取其中第一次出现的染色体来代表具有重复注释的基因。
genes <- genes[!duplicated(genes$ENTREZID),]将基因注释的数据框加入数据对象，数据即被整洁地整理入一个DGEList对象中，它包含原始计数数据和相关的样品信息和基因注释。
x$genes <- genes3、数据预处理
更深的测序总会产生更多的序列。使用edgeR中的cpm函数将原始计数转换为CPM和log-CPM值，来消除测序深度所导致的差异。
lcpm <- cpm(x, log=TRUE, prior.count=2)
在任何样本中都没有足够高的表达量的基因应该从下游分析中过滤掉。edgeR包中的filterByExpr函数提供了自动选取过滤低表达基因阈值的方法。此函数默认选取最小的组内的样本数量为最小样本数，保留至少在这个数量的样本中有10个或更多序列片段计数的基因。
x <- x[keep.exprs,, keep.lib.sizes=FALSE]
下游分析假设所有样品表达值的范围和分布都相似，因此我们需要进行归一化来确保整个实验中每个样本的表达分布都相似。使用edgeR中的calcNormFactors函数，用TMM方法进行归一化（normalization）。当我们使用DGEList对象时，这些归一化系数被自动存储在x$samples$norm.factors。
用limma中的plotMDS函数绘制MDS图，使用无监督聚类来检查样品是否在组内相似而在组间不相似。
par(mfrow=c(1,2))
col.group <- group
levels(col.group) <- brewer.pal(nlevels(col.group), "Set1")
col.group <- as.character(col.group)
col.lane <- lane
levels(col.lane) <- brewer.pal(nlevels(col.lane), "Set2")
col.lane <- as.character(col.lane)
plotMDS(lcpm, labels=group, col=col.group)
title(main="A. Sample groups")
plotMDS(lcpm, labels=lane, col=col.lane, dim=c(3,4))
title(main="B. Sequencing lanes")
在这个数据集中，可以看出样本在维度1和2能很好地按照实验分组聚类，随后在维度3按照测序道（样品批次）分离（如下图所示）。
此外，也可以使用Glimma包的glMDSPlot函数生成互动性MDS图。其左侧面板为一张MDS图，而右侧面板包含一张展示了各个维度所解释的方差比例的柱形图。点击柱形图中的柱可切换MDS图绘制时所使用的维度，且将鼠标悬浮于单个点上可显示相应的样本标签。也可切换配色方案，以突显不同细胞类型或测序泳道（批次）。由下方代码生成的此数据集的交互式MDS图可以在这个网址查看：https://bit.ly/2MIqrbp
groups=x$samples[,c(2,5)], launch=FALSE)
4、差异表达分析
在此研究中，我们想知道哪些基因在我们研究的三组细胞之间以不同水平表达。在我们的分析中，假设基础数据是正态分布的，为其拟合一个线性模型。在此之前，需要创建一个包含细胞类型以及测序泳道（批次）信息的设计矩阵。
design <- model.matrix(~0+group+lane)colnames(design) <- gsub("group", "", colnames(design))
对于一个给定的实验，通常有几种等价的方法可以创建一个合适的设计矩阵。比如说，~0+group+lane去除了第一个因子group的截距，但第二个因子lane的截距被保留。此外也可以使用~group+lane，来自group和lane的截距均被保留。我们在此分析中选取第一种模型，因为在没有group的截距的情况下能更直截了当地设定模型中的对比。
用于细胞群之间成对比较的对比可以在limma中用makeContrasts函数设定。
BasalvsLP = Basal-LP,
BasalvsML = Basal - ML,
LPvsML = LP - ML,
levels = colnames(design))
limma的线性模型方法的核心优势之一便是其可以适应不同的实验复杂程度。无论这篇教程中关于细胞类型和批次的简单实验设计，还是更复杂的因子设计和含有交互作用项的模型，都能够得到相对简单的处理。
对于RNA-seq原始计数或其log-CPM值，方差与均值并不独立，而使用负二项分布来模拟计数的方法假设均值与方差间具有二次的关系。在limma中，假设log-CPM值符合正态分布，并使用由voom函数计算得到的精确权重来调整均值与方差的关系，从而对log-CPM值进行线性建模。
下图左侧展示了这个数据集log-CPM值的均值-方差关系。通常而言，voom图中均值随着方差增加而递减。生物学差异高的实验通常会有更平坦的趋势，其方差值在高表达处稳定。生物学差异低的实验更倾向于急剧下降的趋势。
此外，如果对于低表达基因的过滤不够充分，在图上表达低的一端可以观察到方差水平的下降。如果观察到了这种情况，应当回到最初的过滤步骤并提高用于该数据集的表达阈值。
par(mfrow=c(1,2))v <- voom(x, design, plot=TRUE)
vfit <- lmFit(v, design)
vfit <- contrasts.fit(vfit, contrasts=contr.matrix)
efit <- eBayes(vfit)
plotSA(efit, main="Final model: Mean-variance trend")
limma的线性建模使用lmFit和contrasts.fit函数进行。它们会单独为每个基因的表达值拟合一个模型。然后，通过利用所有基因的信息来进行经验贝叶斯调整，这样可以获得更精确的基因水平的变异程度估计。
为快速查看差异表达水平，以表格列出各对比中显著上调和下调的基因数量。
summary(decideTests(efit))## BasalvsLP BasalvsML LPvsML## Down 4648 4927 3135
## NotSig 7113 7026 10972
## Up 4863 4671 2517
如果需要检测差异表达基因之间的表达量差异倍数是否大于某一阈值（比如本文中的lfc=1，即两倍），可以使用treat方法。
dt <- decideTests(tfit)
summary(dt)## BasalvsLP BasalvsML LPvsML
## Down 1632 1777 224
## NotSig 12976 12790 16210
## Up 2016 2057 190
在多种比较中皆差异表达的基因可以从decideTests的结果中提取，其中的0代表不差异表达的基因，1代表上调的基因，-1代表下调的基因。write.fit函数可用于将三个比较的结果提取并写入到单个输出文件。
vennDiagram(dt[,1:2], circle.col=c("turquoise", "salmon"))
write.fit(tfit, dt, file="results.txt")
使用topTreat函数可以列举出使用treat得到的结果中靠前的DE基因（对于eBayes的结果可以使用topTable函数）。默认情况下，topTreat将基因按照校正p值从小到大排列，并为每个基因给出相关的基因信息、log-FC、平均log-CPM、校正t值、原始及经过多重假设检验校正的p值。
head(basal.vs.lp)## ENTREZID SYMBOL TXCHROM logFC AveExpr t P.Value adj.P.Val
## 12759 12759 Clu chr14 -5.46 8.86 -33.6 1.72e-10 1.71e-06
## 53624 53624 Cldn7 chr11 -5.53 6.30 -32.0 2.58e-10 1.71e-06
## 242505 242505 Rasef chr4 -5.94 5.12 -31.3 3.08e-10 1.71e-06
## 67451 67451 Pkp2 chr16 -5.74 4.42 -29.9 4.58e-10 1.74e-06
## 228543 228543 Rhov chr2 -6.26 5.49 -29.1 5.78e-10 1.74e-06
## 70350 70350 Basp1 chr15 -6.08 5.25 -28.3 7.27e-10 1.74e-06basal.vs.ml <- topTreat(tfit, coef=2, n=Inf)
head(basal.vs.ml## ENTREZID SYMBOL TXCHROM logFC AveExpr t P.Value adj.P.Val
## 242505 242505 Rasef chr4 -6.53 5.12 -35.1 1.23e-10 1.24e-06
## 53624 53624 Cldn7 chr11 -5.50 6.30 -31.7 2.77e-10 1.24e-06
## 12521 12521 Cd82 chr2 -4.69 7.07 -31.4 2.91e-10 1.24e-06
## 20661 20661 Sort1 chr3 -4.93 6.70 -30.7 3.56e-10 1.24e-06
## 71740 71740 Nectin4 chr1 -5.58 5.16 -30.6 3.72e-10 1.24e-06
## 12759 12759 Clu chr14 -4.69 8.86 -28.0 7.69e-10 1.48e-0
5、为结果绘图
plotMD函数可用于绘制均值-差异（MD）图，其中展示了线性模型拟合所得到的log-FC与log-CPM均值间的关系，而差异表达的基因会被重点标出。
xlim=c(-8,13))
Glimma的glMDPlot函数可绘制交互式的均值-差异图，其左侧为与plotMD的输出类似的MD图，右侧为选中基因在各个样本的log-CPM值，下方为结果的表格，可以使用提供的注释（比如基因名）搜索特定基因。
side.main="ENTREZID", counts=lcpm, groups=group, launch=FALSE)
上方指令生成的均值-差异图可以在这里查看：https://bit.ly/2MGbzKF
使用gplots包的heatmap.2函数，我们为basal与LP的对照中前100个DE基因（按调整p值排序）绘制热图。
basal.vs.lp.topgenes <- basal.vs.lp$ENTREZID[1:100]
i <- which(v$genes$ENTREZID %in% basal.vs.lp.topgenes)
mycol <- colorpanel(1000,"blue","white","red")
heatmap.2(lcpm[i,], scale="row",
labRow=v$genes$SYMBOL[i], labCol=group,
col=mycol, trace="none", density.info="none",
margin=c(8,6), lhei=c(2,10), dendrogram="column")
7、使用camera的基因集检验
从http://bioinf.wehi.edu.au/software/MSigDB/以RData对象格式下载Broad Institute的MSigDB c2中的适应小鼠的c2基因表达特征，并运行camera。此外，对于人类和小鼠，来自MSigDB的其他有用的基因集也可从此网站获取。
load(system.file("extdata", "mouse_c2_v5p1.rda", package = "RNAseq123"))idx <- ids2indices(Mm.c2,id=rownames(v))
cam.BasalvsLP <- camera(v,idx,design,contrast=contr.matrix[,1])
cam.BasalvsML <- camera(v,idx,design,contrast=contr.matrix[,2])
cam.LPvsML <- camera(v,idx,design,contrast=contr.matrix[,3])
camera函数通过比较假设检验来评估一个给定基因集中的基因是否相对于不在集内的基因而言在差异表达基因的排序中更靠前。我们选取与我们的数据使用了相同细胞类型的基因集，绘制条码图（barcodeplot）。
index2=idx$LIM_MAMMARY_LUMINAL_MATURE_DN, main="LPvsML")
limma中也有其他的基因集检验，比如mroast的自包含检验。简单来说，camera更适用于在大型数据库中搜寻具有意义的基因集，而mroast测试的是已经确定有意义的一个或少个基因集的显著性。
此外，我们也推荐使用EGSEA进行基因集检验，结合12种优秀的的基因集测试算法的结果来获得生物学结果的整体排名，并且可以生成结果报告网页。如果感兴趣的话，可以搜索EGSEA123查看其使用教程（目前只有英文版，但我们也有在今后将其翻译成中文的打算）。
8、参考文献
Bioconductor Core Team. 2016a. Homo.sapiens: Annotation package for the Homo.sapiens object. https://bioconductor.org/packages/release/data/annotation/html/Homo.sapiens.html.
———. 2016b. Mus.musculus: Annotation package for the Mus.musculus object. https://bioconductor.org/packages/release/data/annotation/html/Mus.musculus.html.
Durinck, S., Y. Moreau, A. Kasprzyk, S. Davis, B. De Moor, A. Brazma, and W. Huber. 2005. “BioMart and Bioconductor: a powerful link between biological databases and microarray data analysis.” Bioinformatics 21:3439–40.
Durinck, S., P. Spellman, E. Birney, and W. Huber. 2009. “Mapping identifiers for the integration of genomic datasets with the R/Bioconductor package biomaRt.” Nature Protocols 4:1184–91.
Huber, W., V. J. Carey, R. Gentleman, S. Anders, M. Carlson, B. S. Carvalho, H. C. Bravo, et al. 2015. “Orchestrating High-Throughput Genomic Analysis with Bioconductor.” Nature Methods 12 (2):115–21. http://www.nature.com/nmeth/journal/v12/n2/full/nmeth.3252.html.
Law, C. W., Y. Chen, W. Shi, and G. K. Smyth. 2014. “Voom: Precision Weights Unlock Linear Model Analysis Tools for RNA-seq Read Counts.” Genome Biology 15:R29.
Liao, Y., G. K. Smyth, and W. Shi. 2013. “The Subread Aligner: Fast, Accurate and Scalable Read Mapping by Seed-and-Vote.” Nucleic Acids Res 41 (10):e108.
———. 2014. “featureCounts: an Efficient General-Purpose Program for Assigning Sequence Reads to Genomic Features.” Bioinformatics 30 (7):923–30.
Lim, E., D. Wu, B. Pal, T. Bouras, M. L. Asselin-Labat, F. Vaillant, H. Yagita, G. J. Lindeman, G. K. Smyth, and J. E. Visvader. 2010. “Transcriptome analyses of mouse and human mammary cell subpopulations reveal multiple conserved genes and pathways.” Breast Cancer Research 12 (2):R21.
Liu, R., K. Chen, N. Jansz, M. E. Blewitt, and M. E. Ritchie. 2016. “Transcriptional Profiling of the Epigenetic Regulator Smchd1.” Genomics Data 7:144–7.
Liu, R., A. Z. Holik, S. Su, N. Jansz, K. Chen, H. S. Leong, M. E. Blewitt, M. L. Asselin-Labat, G. K. Smyth, and M. E. Ritchie. 2015. “Why weight? Combining voom with estimates of sample quality improves power in RNA-seq analyses.” Nucleic Acids Res 43:e97.
McCarthy, D. J., and G. K. Smyth. 2009. “Testing significance relative to a fold-change threshold is a TREAT.” Bioinformatics 25:765–71.
Ritchie, M. E., B. Phipson, D. Wu, Y. Hu, C. W. Law, W. Shi, and G. K. Smyth. 2015. “limma Powers Differential Expression Analyses for RNA-Sequencing and Microarray Studies.” Nucleic Acids Res 43 (7):e47.
Robinson, M. D., D. J. McCarthy, and G. K. Smyth. 2010. “edgeR: A Bioconductor Package for Differential Expression Analysis of Digital Gene Expression Data.” Bioinformatics 26:139–40.
Robinson, M. D., and A. Oshlack. 2010. “A Scaling Normalization Method for Differential Expression Analysis of RNA-seq data.” Genome Biology 11:R25.
Sheridan, J. M., M. E. Ritchie, S. A. Best, K. Jiang, T. J. Beck, F. Vaillant, K. Liu, et al. 2015. “A pooled shRNA screen for regulators of primary mammary stem and progenitor cells identifies roles for Asap1 and Prox1.” BMC Cancer 15 (1). BioMed Central:221.
Smyth, G. K. 2004. “Linear Models and Empirical Bayes Methods for Assessing Differential Expression in Microarray Experiments.” Stat Appl Genet Mol Biol 3 (1):Article 3.
Subramanian, A., P. Tamayo, V. K. Mootha, S. Mukherjee, B. L. Ebert, M. A. Gillette, A. Paulovich, et al. 2005. “Gene Set Enrichment Analysis: A Knowledge-Based Approach for Interpreting Genome-Wide Expression Profiles.” Proc Natl Acad Sci U S A 102 (43):15545–50.
Wu, D., E. Lim, F. Vaillant, M. L. Asselin-Labat, J. E. Visvader, and G. K. Smyth. 2010. “ROAST: rotation gene set tests for complex microarray experiments.” Bioinformatics 26 (17):2176–82.
Wu, D., and G. K. Smyth. 2012. “Camera: a competitive gene set test accounting for inter-gene correlation.” Nucleic Acids Res 40 (17):e133.
更为详尽的信息请点击阅读原文
2019年，遇见更好的自己